% This script conducts market simulations, looking at market efficiency and
% generated behaviors. It simulates a market for T periods, assuming the
% true state is V=1.

clear all;
load('pt_parameters');
load('actual_beliefs');

nonbayesian = 0; % set to use non-Bayesian updating
modalonly = 0; % set to use only modal type
R = 5000; % number of market replications
T = 40; % number of time periods in a market
q = 0.7; % signal quality
quantiles = [0.9,0.5,0.1]; % price quantiles to plot

if modalonly == 1
    [unused,coltouse] = find(params(3,:)==max(params(3,:)));
    params=params(:,coltouse);
end
numtypes = size(params,2);
numtraders = sum(params(3,:));

numup = -4:4; % public signal differences in experiment
pricethreshold_int = q.^numup./(q.^numup+(1-q).^numup); % public beliefs in experiment
% add public beliefs of 0 and 1 to ends for extrapolation
% beyond public beliefs in experiment; in this case, assume
% Bayesian and non-Bayesian agree at 0 and 1
pricethreshold = [0 pricethreshold_int 1];
nonbayesianbeliefs = [0 0; actualbeliefs_nb; 1 1];

parr = zeros(2,R,T+1); % price paths for each set of preferences and replication

for m=1:2 % first pass is risk-neutral, second is for distribution of preferences
    
    % extract parameters for different types
    % in order to break indifference conditions, need to perturb the parameters
    % slightly; for herding make alpha slightly smaller and for contrarian,
    % make delta slightly smaller
    % these perturbations give preferences consistent with the experimental
    % data because the parameters were set at the maximum values consistent
    % with the data
    if m==2
        delta = (params(1,:)>=0)+(params(1,:)<0).*(params(1,:)+1-1E-6); % delta is 1 for herding/risk-neutral types, else delta-alpha+1
        alpha = (params(1,:)<=0)+(params(1,:)>0).*(1-params(1,:)-1E-6); % alpha is 1 for contrarian/risk-neutral types, else 1-(delta-alpha)
        lambda = params(2,:);
    else
        delta = ones(1,numtraders);
        alpha = ones(1,numtraders);
        lambda = ones(1,numtraders);
    end

    for r=1:R % loop over markets
        b = zeros(1,T+1); % public beliefs in each period
        p = zeros(1,T+1); % prices in each period
        b(1) = 0.5;
        p(1) = 0.5;
        signals = rand(1,T)<=q; % random private signals

        for t=1:T % loop over time periods

            % calculate optimal decision for each type and each signal

            % private beliefs for each signal
            if nonbayesian==1 && m==2
                % we have estimated non-Bayesian beliefs for each public belief
                % in the experiment and each private signal - find the closest
                % public belief in the experiment and interpolate

                [unused,closestidx] = min(abs(b(t)-pricethreshold));
                closestbelief = pricethreshold(closestidx);
                if b(t) >= closestbelief
                    loweridx = closestidx;
                    upperidx = closestidx+1;
                else
                    loweridx = closestidx-1;
                    upperidx = closestidx;
                end
                bplus = nonbayesianbeliefs(loweridx,1)+(b(t)-pricethreshold(loweridx))/(pricethreshold(upperidx)-pricethreshold(loweridx))*(nonbayesianbeliefs(upperidx,1)-nonbayesianbeliefs(loweridx,1));
                bminus = nonbayesianbeliefs(loweridx,2)+(b(t)-pricethreshold(loweridx))/(pricethreshold(upperidx)-pricethreshold(loweridx))*(nonbayesianbeliefs(upperidx,2)-nonbayesianbeliefs(loweridx,2));
            else % Bayesian beliefs
                bplus = b(t)*q/(b(t)*q+(1-b(t))*(1-q));
                bminus = b(t)*(1-q)/(b(t)*(1-q)+(1-b(t))*q);
            end
            aplusvec = []; % vectors of actions, one per trader
            aminusvec = [];
            for i = 1:numtypes
                EUbuyplus = EU_la([bplus;1-bplus],[1-p(t);-p(t)],delta(i),alpha(i),lambda(i));
                EUbuyminus = EU_la([bminus;1-bminus],[1-p(t);-p(t)],delta(i),alpha(i),lambda(i));
                EUsellplus = EU_la([bplus;1-bplus],[p(t)-1;p(t)],delta(i),alpha(i),lambda(i));
                EUsellminus = EU_la([bminus;1-bminus],[p(t)-1;p(t)],delta(i),alpha(i),lambda(i));
                % assume trader follows signal if indifferent 
                % action is 1 for buy, 0 for no trade, -1 for sell
                aplus = double((EUbuyplus>=0)&&(EUbuyplus>=EUsellplus)) - double((EUsellplus>0)&&(EUsellplus>EUbuyplus));
                aminus = double((EUbuyminus>0)&&(EUbuyminus>EUsellminus)) - double((EUsellminus>=0)&&(EUsellminus>=EUbuyminus));
                % expand to each trader
                aplusvec = [aplusvec aplus*ones(1,params(3,i))];
                aminusvec = [aminusvec aminus*ones(1,params(3,i))]; 
            end

            % randomly choose a trader and set trade according to realized
            % signal
            randt = randi(numtraders);
            if signals(t)==1
                action = aplusvec(randt);
            else
                action = aminusvec(randt);
            end

            % update beliefs
            % probability of seeing the action conditional on V is sum over
            % types that take action with particular signal times that
            % probability of that signal conditional on V

            % no need to divide through by number of subjects as it cancels out
            % in belief calculation
            praction1 = q*sum(aplusvec==action)+(1-q)*sum(aminusvec==action);
            praction0 = (1-q)*sum(aplusvec==action)+q*sum(aminusvec==action);

            b(t+1) = b(t)*praction1/(b(t)*praction1+(1-b(t))*praction0);

            % update price
            p(t+1) = b(t+1);

        end

        % store market results
        parr(m,r,:) = p;
    end
end

% plot prices
figure;
if modalonly==1
    title('Modal Only');
else
    title('Full Distribution');
end
hold on;
for j=1:size(quantiles,2)
    if j==1
        lstyle = ':b';
    elseif j==2
        lstyle = ':g';
    else
        lstyle = ':r';
    end
    plot(1:(T+1),squeeze(quantile(parr(1,:,:),quantiles(j))),lstyle,'Linewidth',1);
end
for j=1:size(quantiles,2)
    if j==1
        lstyle = '-b';
    elseif j==2
        lstyle = '-g';
    else
        lstyle = '-r';
    end
    plot(1:(T+1),squeeze(quantile(parr(2,:,:),quantiles(j))),lstyle,'Linewidth',2);
end
legend('Benchmark, 9th decile','Benchmark, median','Benchmark, 1st decile','Model, 9th decile','Model, median','Model, 1st decile');
xlabel('Time');
ylabel('Price');
xlim([1 T+1]);

disp('Median pricing error at t=20 (benchmark and model)');
median(1-parr(1,:,21))
median(1-parr(2,:,21))





